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Abstract 

We discuss two methods for exploring energy diffusion in lattices with finite temperature in this paper. 
The first one is the energy-kick (EK) method. To apply this method, one adds an external energy kick 
to a particle in the lattice, and tracks its evolution by evolving the kicked system. The second one is 
the fluctuation-correlation (FC) method. This method is presented recently by one of the present authors 
[Zhao, Phys. Rev. Lett. 86, 11003 (2006)]. In present paper, the formula for calculating the probability 
density function (PDF) using the canonical ensemble is slightly revised and extended to the microcanonical 
ensemble. To apply the FC method, one tracks the motion of the energy initially localized at a small region 
by a properly constructed correlation function of energy fluctuations. Both methods can obtain a PDF of 
energy diffusion. However, we show that the FC method has advantages over the EK method theoretically 
and technically. Theoretically, the PDF obtained by the FC method reveals the diffusion processes of the 
inner energy while the PDF obtained by the EK method represents that of the kick energy. The diffusion 
processes of the inner energy and the external energy added to the system, i.e., the kick energy, may be 
different quantitatively and even qualitatively depending on models. To show these facts, we study not only 
the equilibrium systems but also the stationary nonequilibrium systems. Examples showing that the inner 
energy and the kick energy may have different diffusion behavior are reported in both cases. Technically, 
since applying the energy fluctuations of particles in the system, a set of independent realizations to the 
ensemble average can be achieved by one round of evolution of the system when applying the FC method. 
This advantage enables us to study the long-time diffusion processes in large-scale systems and thus avoids 
the finite-time effect. 
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I. INTRODUCTION 



Diffusion is one of the most important types of motion in nature. According to the time- 
dependent behavior, (r 2 (t)) ~ t a , of the mean square displacement of a conserved quantity, it is 
classified as subdiffusion (a < 1), normal diffusion (a = 1) and superdiffusion (a > 1) [Q,Q]. 
Among various quantities, the diffusion of particle and the diffusion of energy have particular 
importance both for theoretical studies and applications. The diffusion of particles represents 
the transport behavior of mass, while the diffusion of energy determines the transport of heat. 
In studying the particle diffusion, one can directly calculate the PDF by tracking the motion of 
particles. This method has been widely used yi |4|]. However, this method is not available for 
energy diffusion, because a particle can be tagged and tracked while energy may not be. Therefore, 
developing methods to track energy diffusion is highly desirable. This is particularly important for 
a lattice system since particles in the system oscillate around their stationary positions and thus is 
not a diffusion process. 

The key task of characterizing a diffusion process is to obtain the probability density function 
(PDF) of the related quantity, by which one can calculate the time-dependent behavior of the mean 
square displacement and thereby obtain the diffusion exponent a. This exponent is not only served 
to classify diffusion types but also to check certain theoretical predictions. In the past decade, 
the heat conduction problem in low-dimensional systems has attracted intensive attentions [|5|,|6j]. 
It is found that some one-dimensional lattices may show anomalous heat conduction behavior 
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Q7HT1J]- The thermoconductivity k in such a system depends on the lattice size iV as k ~ A r/3 . 
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While some other one-dimensional lattices show normal heat conduction behavior 
recent years, several groups investigated the relation between the exponents a and f3 1117 - 
derived equations to describe their relationship. Establishing deterministic relation between the 
two exponents has theoretical importance. Actually, this attempt is stimulated by the celebrated 
equation obtained by Einstein [12311 who derived the deterministic relation between the diffusion 
coefficient and the coefficient of viscosity, correlating the two different irreversible processes. 
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However, up to now, different relationship equations of a and (3 coexist QJ/7J, 
describes the correct relationship, and even whether the correct one has been found, are still open 
problems. To answer these questions, the exponents a and (3 should be calculated with sufficient 
precision after getting rid of the size effect or finite-time effect. As shown in Ref. [24-26], the 
size effect of (3 disappears in the FPU-/? lattices until N ~ 10 3 — 10 4 . This fact reminds one that 
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the diffusion time should be long enough to avoid the finite-time effect on the exponent a. 

Two different methods have been presented for calculating the PDF of energy diffusion. The 



first one is a straightforward method n 1 9l. \2(X l22l I27M29I1 . The idea is to add a high-energy kick to 



a particle at a fixed position and then, after a period and at other positions, calculate the ensemble 
average of the difference between the current energy density and the average energy density before 
the kick. It is clear that the ensemble average represents the amount of kick energy transported 
to these positions. We call it the energy-kick (EK) method. The EK method indeed manifests the 
diffusion of the kick energy. The second way is presented by one of the present authors, H. Zhao, 
recently [21]. The basic idea is to track the motion of a part of energy by investigating a properly 
constructed correlation function of energy fluctuations. We call it the fluctuation-correlation (FC) 
method. Unlike the EK method, it studies the energy fluctuations of particles in the system. Some 
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groups have realized that this method is more reasonable than the EK method II22L 

The purpose of this paper is, on one hand, to make a detailed comparison between the EK 
method and the FC method. We argue that the two methods indeed explore different diffusion 
processes. That is to say, the EK method describes the relaxation process of the kick energy or the 
diffusion behavior of the energy externally added to the system, while the FC method represents 
the diffusion behavior of the inner energy or the intrinsic diffusion behavior of the energy in the 
system. In certain systems, the memory of the history will be lost quickly with time evolution. 
In this case, the diffusion of the kick energy and the inner energy has no qualitative difference, 
though the PDFs obtained by the two method still have quantitative difference since the diffusion 
of the kick energy depends on the amplitude of kicks. In other systems, the memory of the his- 
tory is not totally lost. In this case, the kick energy and the inner energy may include different 
information. For most lattices in equilibrium with uniform temperature, the kick energy and the 
inner energy diffuse in a symmetric way, and the PDFs obtained by both methods are qualitatively 
similar. These lattices usually have symmetric interaction potentials. While for lattices with asym- 
metric interaction potentials, the inner energy and the kick energy may diffuse in different ways. 
The PDFs obtained by the FC method are symmetric and those obtained by the EK method may 
be asymmetric. In nonequilibrium lattices with stationary temperature gradients, even for those 
lattices with symmetric interaction potentials, the PDFs obtained by the two methods may be dif- 
ferent qualitatively. Therefore, we consider not only the energy diffusion in equilibrium systems 
but also in nonequilibrium systems. 

On the other hand, we further explain the rationale of the FC method and describe the technical 
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details for applying this method. We present formulae for the canonical ensemble and the micro- 
canonical ensemble respectively to calculate PDFs of energy diffusion. A formula suitable for the 
canonical ensemble has been presented in [2l|. The one presented here is slightly different from 
the old one. We explain that the new formula is more reasonable. Also we emphasize that the FC 
method has great technique advantage over the EK method. The technical disadvantage of the EK 
method has been pointed out in Ill9n . The FC method can shorten the computation time by several 
orders comparing with that with the EK method. 

The paper is managed as follows. In next section we introduce the two methods for studying 
energy diffusion in equilibrium one-dimensional lattices. The rationale of the FC method is ex- 
plained and the strategy for applying this method is described. An example showing the qualitative 
difference between the inner energy and the kick energy is reported. Section III is contributed to 
investigating the energy diffusion in stationary nonequilibrium lattices, i.e., in one-dimensional 
lattices with stationary heat flux and stationary temperature gradients. This section further reveals 
the diffusion behavior difference between the inner energy and the kick energy. The last section is 
the conclusion and the discussion of the paper. 



H. ENERGY DIFFUSION IN EQUILIBRIUM LATTICES 

A one-dimensional lattice is usually described by the Hamiltonian 

N 2 

H = Y^Hi, H = ^ + U{ Xi ) + V[x % - x f _i), 

i=l 

where U (xj) denotes the on-site potential and V(x i+ i — Xi) the interaction potential. In the present 
paper we study two types of the lattice models, without or with on-site potentials. The illustrating 
examples include the FPU-/3 model with 



p 2 1 1 
Hi = H — (x{ — Xj_i) 2 H — (xi — Xi-i^ 4 
2m 2 4 V 



and the lattice (b A model with 



p 2 I 1 



2 , 

2m ' 2 V ~' ~ L ~ 1J 4* 



Without loss of generality, we set mass m = 1, and the Boltzmann constant = 1 in following 
discussions. Fixed boundary condition is applied for both models. We also employ occasionally 
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FIG. 1: (color online). The PDFs of energy diffusion at T = 0.5. (a) The PDFs at t = 50 obtained by the 
EK method (A) and the improved EK method (B) for the FPU-/? model, (b) The PDFs at t = 100 obtained 
by the FC method(A), and the EK method with A^(0) = 5E (B) and A#;(0) = 10E (C) for the FPU-/? 
model, (c) The PDFs at t = 200 for the lattice </> 4 model. A, B and C are same as in (b). (d) The PDFs 
at t = 100 for the FPU-/3 model obtained by the FC method using microcanonical ensembles with eqI2^4) 
and with eq|4](i?). The line C gives p(r, t) = —1/N for reference. 

2 2 

the harmonic model with H{ = P- + \{xi — Xj_i) 2 , the Toda model with Hi = j*- + e~^ Xi ~ Xi ~ 1 ' + 
(xi — — 1 and the quartic-FPU model with if, = ^- + Mxi — Xj_i) 4 for special purposes. 

When adding heat baths to the two ends of a lattice one obtains a canonical system. A canonical 
system is an open system and thus the total energy of the system is not conservative. When heat 
the lattice to a finite temperature and then take away the heat bath, one obtains a microcanonical 
system. A microcanonical system is an isolated system with total energy conserved. In both 
cases, after the stationary state is achieved, every particle has identical temperature, T = (pf), and 
identical average energy, E = (Hi). 

The problem to be investigated is: how does the energy Hi or the energy difference AHi = 
Hi — E initially located at the ith particle spread out over the lattice as a function of time? To 
find the solution to this problem, the EK method adds a kick Af7j(0) = Hi(0) — i/j(0) to the ith 
particle at t — 0, and calculate the part AHj(t) transported to the jth particle at time t, where 
Hi(0) is the energy before the kick while Hi(0) represents the energy after the kick. It is clear that 
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without the external kick the ensemble average (^Hj(t)j should be equal to the average energy E. 
Therefore, the ensemble average (AHj(t)) = /Hj(t)\ — E can be explained as the average part 
of the kick energy transported to the jth particle at time t. One can thereby define 

(AH jit)) 

p(r,t) = ^- z f (1) 

A#i(0) 



as the PDF to describe the probability of the kick energy transported to the position r at time t, 
where r = j — i. 

As an example, in Fig. [TJa) we show the PDF obtained by the EK method for the FPU-/3 
model with temperature T = 0.5 and at t = 50. The ensemble average is over 4.2 x 10 5 different 
realizations. For each realization the lattice is developed to t = 50, starting by adding an external 
kick to the middle particle. For the sake of simplicity we fix the amplitude of kicks at Aifj(0) = 
5E. The figure indicates that the profile of the PDF is already distinguishable, but the fluctuations 
are still large. 

A different strategy can be applied to improve the convergence of the PDF for the EK method 



[221. The idea is to copy the system before the kick as a reference system. Then add the kick 
to the middle particle of the original system and evolve the pair of the systems simultaneously to 
compute the energy difference between the two systems by AHj(t) = Hj(t) — Hj(t), where the 
first term in the r.h.s is the energy of the jth particle of the kicked system while the second is that of 
the reference system. It is clear that the difference represents exactly the kick energy transported 
to the jth particle. By averaging the same amount of realizations as in above we obtained the PDF 
and show it also in Fig. Etta). One can find that the fluctuations are dramatically smoothed. The 
PDFs obtained by the EK method hereafter are calculated using this strategy. 

The EK method is presented in the spirit of the linear response theory[31], which demands the 
kick energy small enough to keep a linear response. However, practically, it can not be too small 
otherwise the fluctuations of the ensemble average will be too large. In this situation, the PDF will 
depend on the amplitude of the kick energy. In Fig. [jjb) we show the PDFs calculated with kicks 
of amplitude Aifj(0) = 5E and of amplitude Aifj(0) = 10E. It is obvious that the kick energy 
with AHi(0) = 10E diffuses faster than that with AHi(0) = 5E. The same effect remains in the 
lattice 4 model, as Fig. djc) shows. Nevertheless, in this model the high-energy kick diffuses 
slower than that of the low-energy kick, which confirms the observation that the higher of the 



energy, the stronger of the localization effect in this model 132(1 . Therefore, the EK method studies 
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FIG. 2: (color online). The initial (spatial) correlation (AHi(0)AHj(0)} of the FPU-/3 model (a) and the 
lattice (f) 4 model (b). 

indeed the diffusion process of the kick energy; it is not directly explore the intrinsic diffusion 
behavior in the system. 

The desirable way is to study directly the inner energy, i.e., the energy initially localized at 
a particle or in a small region in the lattice. The difficult is that one can not directly track this 
part of energy because it can not be marked, as does in the case of particle diffusion. Let AfZj(0) 
represents the energy fluctuation initially located at the ith particle. Then if a portion of Aifj(O) 
is transported to the jth particle at time t, the energy fluctuation AHj(t) should have correla- 
tion with AfZj(O) physically. The idea of the FC method is to track the motion of the inner 
energy by calculating the correlation of the energy fluctuations. The first step of applying the FC 
method is to determine the part of energy, denoted by AH(0), to be studied. We show how to 
do this by examples. In Fig. Ha) and Ob), we plot (A^(O)A^(O)) of the FPU-/3 model and 
the lattice 4 model respectively, where r = j — i. Canonical systems with N = 501 are em- 
ployed for both models. It can be seen that (AiJ i (0)Ai7 J -(0)) = (within numerical errors) for 
j ^ i and (A^(O)A^(O)) ^ for j = i for the FPU-/3 model, which indicate that A#;(0) 
has no initial correlation with the energy fluctuations of other particles. In this case, we employ 
AH(0) = AHi(0) as the energy to be studied. While in the case of the lattice <p A model, it has 
(AHi(0)AHj(0)) > for several nearby particles, which implies that Ai7j(0) is intrinsically 
correlated with a part of energy besides AHi(0) itself. Because this part of energy can not be 
separated from AH^O), we have to consider it as a whole as the objective AH(0) to be studied. 
We call AH(0) as the adherence energy of AHi(0). 

It is easy to understand why in certain systems Ai^(0) adheres a small part of energy while 
in other systems it does not. In a lattice with finite temperature, at the same moment the fluc- 
tuations of the kinetic energy of particles are independent because the velocity of a particle at a 
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moment puts no restriction on that of other particles. For a class of lattices with only interaction 
potentials, such as the FPU-/3 model, the potential is determined by the relative displacements. A 
configuration of x,- L — which determines the potential part of H, L is completely independent of 
the configuration of x i+ i — %i which determines the potential part of H i+1 . As a result, Hi and 
H i+ i are independent with each other and thus (Aifj(O)Aifj(O)) = for j ^ i. For another class 
of lattices with on-site potentials, such as the lattice 4 model, the situation is different. When 
the ith particle moves to the left/right to its equilibrium position, the neighbor particles will tends 
to follow it because of the interaction potential between them. This effect results in a positive 
correlation, i.e., the bigger the on-site potential xj/A, the bigger the xf ±1 /A, and vice versa. 
We then define 

(Ag.(O)AgjM) 
P[ ' ' (Aft(0)A7/(0)) W 

as the PDF to describe the motion of AH(0). The reason is as follows. First, AH(0) never disap- 
pears before it spreads out of the lattice because of the energy conservation, i.e., it is a conserved 
quantity. Thus, one has £\ A#;(0)A#j(i) = A#;(0)A#(0), which gives f p(r,t)dx = 1. 
Second, AH^O) and AH(0) are positively correlated, and Aifj(t) and AH(0) should also be 
positively correlated since the former is a part of the latter, one obtains p(r, t) ^ 0. Non-negativity 
and normalization conditions are basic requirements to be a PDF. Finally, it is clear that p(r, t) is 
proportional to AH'-(t) / AH(0). Therefore, p(r, t) can be considered as the probability of finding 
AH(0) at the position r at time t, or equally, it represents the rate of AH(0) to be transported to 
r at time t. 

The problem is how to calculate (AiJj(O)AiJj(t)). At time t we divide the energy fluctuation 
at the jth particle into two parts, i.e., AHj(t) = AH'-(t) + AH"(t), where AH'-(t) represents the 
part of AH(0) being transported to the jth particle, and AH"(t) comes from other sources. For 
a canonical system, AH"{t) has no correlation with AiJj(O) since it comes from {other parts of 
the system or} the heat baths. In this case, it should has (AiJ i (0)AiJ ? ' / (t)) = 0. As a result, we 
can calculate (Aif(O)Aifj(i)) by (AH^AH'^t)) = (AH(0)AHj(t)) and (AHi(0)AH(0)) by 
(Aifj(O)AiJ(O)) = J2j {AHi(0)AHj(0)). In other word, in studying canonical systems one can 
apply 

{AimAHM 
P[ ' ' (AHi(0)AH(0)) { ) 

to calculate the PDF of AH(0). 

For a microcanonical system, however, the condition (AHi(0)AH"(t)) = fails because of 



8 



the conservation of the total energy of the system. The conservation of the total energy implies 
£V AHj(t) = 0, which can be rewritten as £\[A#<(t) + AH'!(t)\ = AH(0) + £\ = 0. 

Multiplying AH^O) to the equation one obtains £\ Aifj(0)A.£fj'(0) = -AiJj(0)Aif(0). By 
calculating the ensemble average it appears as £\ (A#;(0) A#j'(£)) = - ( Ai^(0) AH (0)). Sup- 
pose that the lattice includes N particles. Because the nonvanishing correlation (AHi(0)AH"(t)} 
is resulted from the conservation of the total energy, it is an intrinsic correlation of the system 
and should be independent of time. It is reasonable to assume that (Aifj(O) AH" (i)) should have 
identical value for each particle. One then derives 

(AH^AH'jit)) = -i (A^(O)Atf (0)> . 

Notice that AH'Jt) = AHj(t) — AH'Ht), the PDF of the energy diffusion in a microcanonical 
system appears as 

(A^(O)A^.(t)) J_ 
A ' ' (AHi(0)AH(0)) N' K ) 

These two equations[3]and|4]trustily represent the definition equation|2] since one can check that 
the p(r,t) obtained by them satisfy the conditions J p(r,t)dx = 1 and p(r,t) ^ 0. In Fig. [TJb) 
andQJc), we show the p(r, t) calculated by the FC method using canonical ensembles. One can 
see that for either the FPU— (3 model or the lattice </> 4 model the condition p(r, t) ^ is satisfied 
within numerical errors. In Fig. [3l we show J p(r, t)dx as a function of time for the two models 
at T = 0.5, which indicates that J p(r, t)dx = 1 is also satisfied with high precision. 

To obtain an isolated system with given temperature, we heat the model to the stationary state 
with the expected temperature and then take away the heat bath. Applying such a microcanonical 
system to calculate the PDF, one should employ the eq. |4] instead of the eq. [3l In Fig. [Del), the 
dotted line shows the p(r, t) calculated by eq. |4] while the solid line shows that of calculated by 
eq. [3] in the case of the isolated FPU— model. It can be seen that the condition p(r, t) ^ is 
well satisfied in the former case while fails in the latter case. The condition J p(r, t)dx = 1 is also 
satisfied within numerical errors as shown in Fig. [3] in the former case, while it can be checked 
that J p(r, t)dx = in the latter case. Moreover, we have checked that the PDFs obtained with 
the canonical ensemble and the microcanonical ensemble are identical with each other at the same 
temperature within numerical errors. 

From Fig. \V[b) and \V[c) one can realize that a PDF obtained by the EK method will approach 
the corresponding PDF obtained by the FC method when the kick energy Ai^(0) is small enough. 
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FIG. 3: (color online). J p(r, t)dr versus time t for the FPU-/3 model obtained with canonical ensemble(A) 
and microcanonical ensemble (B), and the lattice </> 4 model with canonical ensemble(C) and microcanonical 
ensemble (D) 



However, the smaller the kick, the slower the convergence of the PDF. As can be seen in Fig. 
QIb), with the same amount of the realizations, the fluctuations of the PDF obtained with the kicks 
AHi(0) = 5E is much larger than that of the PDF obtained with A-fiTj(O) = 10E. Therefore, 
decreasing Af/j(0) needs to dramatically increase the amount of the realizations to achieve a 
reliable ensemble average. 

We now describe the technical details in applying the two methods. Let us take the FPU-/5 
model as an example. For this model, as shown in 112 ill , the two soliton-like packets on the PDF 
move at a constant velocity, which is supersonic and depends on the temperature of the system. 
With the dimensionless unit, the soliton-like packets move with v > 1. In this case, supposing that 
the diffusion starts at the middle particle of the lattice and one wants to study the diffusion process 
up to a time length t = t c , the lattice must have a size N > 2t c to keep the energy diffusing within 
the lattice during this period. 

To apply the EK method, one may employ a lattice with length N + 1 and evolve the system 
with a sufficiently long time to relax it to a stationary state. Then copy a reference system before 
adding a kick Aifj(0) to the middle particle of the lattice. And evolve both the reference system 
and the kicked system for a time t c to obtain a set of realization AHj (t) . This step, i.e., kicking and 
evolving the two systems, is repeated again and again to obtain the ensemble average (^AHj(t)^. 
Each round of repeat contributes one realization to the ensemble. In this way, evolving a couple 
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FIG. 4: (color online). Illustration of the strategy of the FC method for calculating the ensemble average. 
Each of the black particles is applied as independent source where the diffusion is starting at. 

of systems for a period t one obtains t/t c realizations to the ensemble. 

To apply the FC method, we apply a lattice with a size N/2 + M + N/2, and select M/£ T 
particles in the middle segment as the sources that the diffusion starts, as illustrated in Fig. SI Here 
£ T is a constant which should be sufficiently bigger than the spatial correlation length between 
particles. For the FPU-/3 model, one can set £ T = 1 since Fig. Eta) has indicated that there is 
no spatial correlations between different particles. For the lattice (p 4 model, Fig. [2tc) has shown 
that the correlation of a particle with its next nearest neighbors has closely approached zero and 
one can set £ T = 3. In this way, each selected particle can be considered as an independent 
source and contributes realizations to the ensemble average independently. Based on the same 
principle, for each of the selected particles, one can further apply a set of continuous records, 
AiZj(O), AHi(t T ), AHi(2t T ), AHi(3t T ),- ■ ■ , as the independent sources, where t T is a constant 
which should be sufficiently bigger than the autocorrelation length of the particle. For the FPU-/3 
and the 4 models, it can be easily checked that one can apply t T = 2 and t T = 3 respectively. In 
this way, each solid circle in Fig. |4]is applied as an independent source of energy diffusion, and 
contributes a realization to the ensemble. 

In this way, by evolving the lattice with a time t one can obtain M(t — t c )/£ T t T realizations to 
the ensemble, which is about Mt c /^ T t T times of those obtained in the case of the EK method if 
t > t c . 

In the FPU-/3 model with temperature T = 0.5, for example, the soliton-like packets on the 
PDF move with a velocity v ~ 1.3. To investigate the diffusion process extended to t c = 1000, 
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one can employ a lattice with iV = 3001 when applying the EK method and apply the middle 
particle as the source that the energy starts to diffuse. While applying the FC method, one may 
employ a lattice with 4000 particles and selects the source particles from the middle segment with 
M = 1000. In both cases, the energy diffusion remains within the lattices during the period t c , 
since within this time the soliton-like packets can spread over about 1300 particles. Evolving the 
two lattices with the same time t, the realizations to the ensemble obtained by the FC method is 
about Mt c /^ T t T ~ 5 x 10 5 times of the realizations obtained by the EK method. 

Figure \5\ (a) shows the PDFs of the FPU-/? model at t c = 800 calculated by averages over 
3 x 10 9 and 3 x 10 11 realizations respectively. It can be seen that the fluctuations of the PDF are 
suppressed to a reasonable level in the latter case while is still remarkable in the former case. Even 
using the FC method, to obtain 3 x 10 11 realizations for such a long time diffusion processes is still 
a hard task for serious computation, to say nothing of the EK method. Fortunately, we achieved 
the goal by parallel computations using 28 CPUs with half a month. If one wants to obtain the 
same amount of realizations by the EK method, he must use about 5 x 10 6 CPUs with the same 
computation time. 

Achieving a sufficiently long diffusion time is necessary to avoid the finite-time effect. The 
finite-time effect can be ignored for the lattice 4 model, in which the exponent a tends to a = 1 
within t c < 100. For other types of lattices with anomalous energy diffusion, such as the FPU-/3 
model and the quartic-FPU model, the finite-time effect is remarkable. The quartic-FPU model 
represents the unharmonic (or high-temperature) limit of the FPU-/3 model. As pointed in [33], 
the qualitative statistical property of this model is temperature-independent because of the scaling 
behavior. In other words, systems with different temperatures can be scaled together by a proper 
scaling transformation. 

As shown in Fig. (2 a), the PDFs of the energy diffusion for the two models are qualitatively 
similar with each other, though the soliton-like packets on the PDF of the quartic-FPU model are 
obviously smaller than that of the FPU-/? model. In Fig. Ob), we plot the (r 2 (t)) as a function of 
time with the log-log scale for the two models. Roughly, it seems that both models show a power- 
law relationship between (r 2 (t)) and t. To explore the finite-time effect, we plot a by sectional 
fitting, as Fig. 6(c) shown. In the plot, each point of a is obtained by fitting 4 data points, i.e., it 
represents the result in a time interval At = 200 correspondingly. One can see that a changes with 
time initially and converges gradually to a constant, and the finite-time-dependent behavior of a 
for the two models have different features. It can be seen also that a approaches 1.41 at t c > 600 
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FIG. 5: (color online). The finite-time effect of the diffusion exponent a. (a) The PDFs at t = 800: A, for 
the quartic-FPU model averaged over 6 x 10 11 realizations; B and C, for the FPU-/3 model averaged over 
3 x 10 9 and and 3 x 10 11 realizations respectively, (b) The log-log scale plot of (r 2 (t)) against t for the two 
models. (c)The time-dependent effect of a for the two models. In both (b) and (c), the squares represent the 
results of the FPU-/3 model and the triangles are for the quartic-FPU model. 

for both models. This result tends to support the relationship equation of a and j3 presented in 
Ref. [18]. However, we still can not sure whether the diffusion exponents have converged exactly 
at t c = 800. To check it one needs to extend the computation to a longer diffusion time, which 
already exceeds the capability of our computer system. At least, this fact reminds one to be careful 
in checking the relationship between a and j3 by numerical calculations. Instead of studying the 
diffusion process with the time scale of t c ~ 100 as done by some previous researchers, a diffusion 
time t c > 700 at least is needed to avoid the finite-time effect. 

According to the above results, it seems that the PDFs obtained by the two methods are al- 
ways symmetric and have only quantitative differences. The symmetric PDFs obtained by the FC 
method represent the intrinsic behavior of the inner energy. Any part of energy must moves in 
symmetric way statistically, otherwise the equilibrium can not be maintained. However, the sym- 
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FIG. 6: (color online). The PDFs of the Toda model at T = 0.5 and t = 100: A, by the EK method with 
positive kicks; B, by the EK method using the same number of positive and negative kicks; C by the FC 
method. 

metric PDFs obtained by the EK method are particular characteristic of a class of special lattices. 
By reconsidering the models discussed above, one can find that they have a common feature, i.e., 
the interaction potential in these models are symmetric. As a consequence, the kick energy added 
on a particle, either with positive or negative momentum, will be equally transported to the op- 
posite directions. In these lattices, the PDFs obtained by the EK method are always symmetric. 
On the contrary, for lattices with asymmetric interaction potentials, the PDF obtained by the EK 
method depends on the way of the kicks. A typical model with asymmetric interaction potential is 
the Toda model. Figure [6] shows the PDFs computed by the FC and the EK methods respectively 
for this model in equilibrium with T = 0.5. Physically, the inner energy must have equal proba- 
bility to move to both sides to maintain the equilibrium, the PDF should be symmetric. The figure 
indicates that the FC method explores this feature correctly. To applying the EK method, we add 
the kicks in two ways. One is to always add kicks with positive momenta. Notice the asymmetric 
feature of the interaction potential in this model, a kick with positive momentum will transport 
more energy to the r.h.s than to the l.h.s. As a result, the PDF of the kick energy should be asym- 
metric, which is confirmed by the figure. Another way is to add positive and negative kicks with 
equal probability. In this case, one can obtain a symmetric PDF as the figure shows. Thus, the 
PDF of the kick energy depends on the way of the kicks. Furthermore, the figure explores another 
serious problem of the EK method. While the condition p(r, t) ^ well-satisfied for the PDF 
obtained by the FC method, it fails for the EK method in certain intervals, and thus fails to be a 
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FIG. 7: (color online). The PDFs at t = 80 calculated in the case of stationary noequilibrium systems by 
the EK method (A) and the FC method (B). (a) The harmonic model, (b) The FPU-/3 model, (c) The lattice 
c/> 4 model. In each model, the lattice size is N = 501 and the temperatures of the heat baths are fixed at 
T+ = 1 and T_ = 0.5. 



probabilistic density function. Therefore, for this model, the EK method is indeed unsuitable for 
studying the diffusion process. 



HI. ENERGY DIFFUSION IN NONEQUILIBRIUM LATTICES 

Studying the diffusion in stationary nonequilibrium lattices provides us more clear examples 
to reveal the different diffusion behavior of the inner energy and the kick energy. When coupled 
with two heat baths with different temperatures, a lattice will approach a nonequilibrium stationary 
state after a relaxation process. A constant heat flux is then established along the lattice. In such a 
nonequilibrium lattice, the temperature as well as the average energy can only be defined locally, 
i.e., Tj = (p?), and E { = (Hi). The AiJj(O) = Hi(0) — E { depends also on the position of the 
particles, and, therefore, the energy diffusion starting at different particles may not be statistically 
identical. In this situation, we study the diffusion starting at the middle particle of a lattice as 
the illustrating example. Consequently, in a nonequilibrium lattice, one can no longer employ a 
segment with M particles to calculate the ensemble average, as doing in the case of equilibrium 
systems. However, the FC method still has the technique advantage in calculating the ensemble 
average since along the time axis the strategy described in Fig. |4]is still available. 

In Fig. [7J we show the PDFs obtained by the EK method (dotted lines) and the FC method 
(solid lines) for the harmonic model, the FPU-/3 model and the lattice </> 4 model. The temperature 
of the heat bath in the l.h.s is T + = 1 and is T_ = 0.5 in the r.h.s. Each lattice includes N = 501 
particles and the diffusion time is t = 80. It can be seen that except of the </> 4 model, the PDFs 
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obtained by the EK method are symmetric while those by the FC method are asymmetric. In more 
detail, it can be checked that a ~ 1 for the PDFs obtained by the EK method. For the PDFs 
obtained by the FC method, it has a ~ 1.4 for the harmonic model, a ~ 1.1 for the FPU-/? model, 
and (7 — 1 for the lattice 4 model. 

The PDFs with o ~ 1 obtained by the EK method represent the diffusion behavior of the kick 
energy. As has been pointed out in the end of last section, for lattices with symmetric interaction 
potentials the kick energy is transported to the opposite directions equally. Either in equilibrium 
or in nonequilibrium, the environment around the middle particle is similar locally, and thus the 
kick energy is divided into two equal parts for this type of lattice. The two parts of energy are 
transported towards the opposite directions hereafter. 

The FC method correctly explores the diffusion behavior of the inner energy. The harmonic 
model is the only one been solved rigidly for the heat conduction problem [7]. The mechanism of 
the heat transport in this model is clear, i.e., the energy delivered by one heat bath will be trans- 
ported to the opposite side directly since there is no interaction among normal modes (phonons) in 
this model. As a result, the energy fluctuation AiJj(O) at the middle particle can be exactly divided 
into two parts Ai^(O) = A/f +(0) + Aifr(O). The part A/f +(0) comes from the l.h.s heat bath 
and will move to the r.h.s, and the part AH[(0) comes from the r.h.s heat bath and will move to 
the l.h.s. In the case of T + = TL, it has AH^(0) = AH~(0) and in a later time the amount of 
the energy Aifj(O) transported to both sides of the middle particle will be the same. While in the 
case of T + > TL, it has AH^(0) > AH^(0) and in a later time the energy moves to the left will 
be more than that moves to the right. This is the expected diffusion property of the inner energy in 
this lattice. We define a parameter a as 

f o °° p(r,t)dx 



a 



to measure the rate of the energy spreading to the two sides of the middle particle. When applying 
the FC method, it is easy to obtain a = ( AHl ^ AH i @>) Q ne jjj US d er i V es cr = 1 in equilib- 

J (Affi(O)Affr(O)) 

rium lattices and a > 1 in nonequilibrium lattices, since AH^(0) = AH[(0) in the former and 
Aif +(0) > A#r(0) in the latter cases. 

In the case of the FPU-/3 model, as pointed out in I21L the soli ton-like packets keep the initial 
memory of the direction wherever they move. Equivalently, in this model the energy in the soliton- 
like packets keeps its initial memory of moving direction. As a result, the energy fluctuation 
Aifj(O) at the middle particle can be divided into three parts, AH^O) = AH+(0) + Aifr(0) + 
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AH®(0). Here AH+(0) and AH~(0) come from the high- and low-temperature heat baths and 
still keep the memory of their moving directions; they will travel to the opposite sides. The part 
AH®(0) represents the no-memory part which has equal probability to move to either direction. 
Because of A//,- > AH- , one can still predict a = -) — — ) ^ — (- > 1. 

1 V 7 * V 7 F (Aif i (0)AHr(o))+(Aif i (0)Aif?(0)) 

For the lattice 4 model, the Gaussian PDF as well as the derived property of (r 2 (t)) ~ t° 
indicates that the energy diffuses in this system normally. Normal diffusion implies a no-memory 
random walk. The energy starting at any point will loss totally the initial memory of its moving 
direction. Thus, Aif;(0) at the middle particle has no memory of the heat baths and will diffuse to 
either side with equal probability. In this case, one should expect a = 1 for this model. 



IV. SUMMARY AND DISCUSSION 



In summary, the FC method reveals the intrinsic diffusion behavior of the energy in the system, 
while the EK method displays the diffusion behavior of the kick energy. The inner energy and the 
kick energy may carry different information, and they may diffuse in different ways. 

In lattices with normal energy diffusion, such as in the lattice 4 model, the energy will totally 
lose the initial memory. In these lattices, the behavior of the kick and the inner energy appear qual- 
itatively similar because both of them will lose the initial information totally. The PDFs obtained 
by the two methods appear as Gaussian functions, either in equilibrium or in nonequilibrium. 
However, the results of the two methods are quantitatively different. The PDFs obtained by the 
EK method depend on the amplitudes of the kick energy, only when the kicks are weak enough 
can it approach the results of the FC method. 

In lattices with anomalous energy diffusion, the energy may keep part of the initial memory. 
In equilibrium case, the memory informtion about the moving directions included in the inner en- 
ergy is symmetric, thus the PDFs obtained by the FC method exhibits such symmetry. When the 
interaction potentials in lattices, such as the harmonic model or the FPU-/3 model, are symmetric, 
the kick energy is thus divided into two equal parts, one moves to the right and the other to the 
left. As a result, the PDFs obtained by the EK method are also symmetric, although they may 
be quantatively different from those obtained by the FC method. When the interaction potential 
is asymmetric, such as in the Toda model, however, the kick energy may diffuse asymmetrically, 
depending on the way how the kicks are added on. Because of the asymmetric potential, a kick 
to the right, for instance, will transport more energy to the r.h.s than to the l.h.s., and vice versa. 
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In this kind of lattices, if one adds positive and negative kicks to the particle with equal prob- 
ability, the PDFs calculated by the EK method are symmetric, otherwise the PDFs must appear 
asymmetrically. 

In nonequilibrium case, the heat baths added on the two ends of a lattice have different tem- 
peratures. As a result, the parts of energy keeping the memory of the corresponding heat bathes 
arrived at the middle particle are different. The PDFs obtained by the FC method are asymmetric 
and reveal this intrinsic feature of inner enrgy, while the EK method fails to do so because the kick 
energy keeps no information of memory of the heat bathes. 

Besides the solid theoretical basis, the FC method has remarkable technique advantage over 
the EK method. Because the FC method uses the energy fluctuations of particles in the system, 
and any independent energy fluctuation can be treated as AiJj(O) to calculate the ensemble aver- 
age, thus one can obtain a large set of independent realizations to the ensemble by one round of 
evolution of the system, as illustrated in Fig. @J This strategy can reduce the computation time 
dramatically. The adventage becomes quite remarkable when one studies the long-time diffusion 
behavior. Comparing with the EK method, the time to evolve the system is shortened by at least 
five orders when one studies a diffusion process extended to t c = 1000. With this strategy, we cal- 
culated the diffusion exponents a of the FPU-/? model and the quartic-FPU model upto t = 800, 
and observed that the size-effect on a obviously exist till at least t = 600. This fact indicates that 
the a calculated in some previous works by the EK method with the diffusion time upto t ~ 100 
is questionable. 

To apply the FC method, we introduced the formulae for calculating the PDF of energy dif- 
fusion available for canonical and micro-canonical ensembles respectively. In the formulae, the 
conserved quantity being investigated is AH(0) instead of Ai/j(0). This is because that in certain 
models the energy AHi(0) adheres constantly a part of energy. To guarantee the normalization 
of p(r,t), we have to consider Aif(0) as the energy to be studied. Therefore, the first step to 
apply the FC method is to detect the energy AH(0) adhered to AHi(0) by checking the correla- 
tion (AHi(0)AHj(0)). The energy AH(0) needs not to be calculated explicitly, but we demand 
it locating in a small region and positively correlated to Aifj(0). When these conditions are ful- 
filled, AH (0) is available for probabilistic descriptions and the p(r, t) describes the probability of 
finding AH(0) at position r and time t. 

The definition of the PDF for canonical system presented in this paper is slightly different from 
that presented in the Ref. 112 ill , in which AfZj(0) is directly applied as the conserved quantity to be 
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studied. This treatment results no problem for certain lattices such as the FPU-/3 model. In such a 
lattice, it has AH(0) = AH^O) because of (AfZi(O)AfZj-(O)) = for i ^ j. For other systems, 
such as the lattice 4 model, AH(0) is slightly bigger than Aifj(O). As a result, the p(r, t) obtained 
can not be normalized exactly. With higher precision, we have checked that J p(r,t)dx ~ 1.07 for 
the lattice 4 model at T = 0.5 when applying the previous definition. Because J p(r, t)dx is still 
time-independent, it will not alter the value of the exponent a. However, as a PDF, normalization 
is an expected feature. It is thus better to consider AH(0) as the conservated quantity to be studied. 
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